
#include "simulator.hpp"
#include "settings.hpp"
#include <cmath>
#include <iostream>

#ifdef CAVITY_MULTITHREADING_AVAILABLE
    #ifndef NDEBUG
    boost::mutex debug_mutex;
    #endif
#endif

simulator::simulator(thread_method & tmeth_, mode_method & mmeth_, grid & g_, grid & lap_,
                     grid::region & reg_, grid::boundary & bound_, uword n_, grid::data_type dt_,
                     uword sub_period_, uword norm_period_)
                : tmeth(tmeth_), mmeth(mmeth_), reg(reg_), bound(bound_), g(g_),
                  lapf(lap_), relaxf(lap_,dt_), n(n_), subtraction_period(sub_period_),
                  normalization_period(norm_period_)
{
    //empty
}

void simulator::operator()() //Here is the simulation
{
    for(uword i = 0; i < n; ++i)
    {
        g.iterate_region(reg,lapf);
        tmeth.wait();	//multithread
        g.iterate_region(reg,relaxf);
        if(i % subtraction_period == 0)	//high mode
        {
            mmeth.subtract(reg);	//high mode
            if(i % normalization_period == 0)
                normalize(1.0);
        }
        mmeth.adjust_boundary(bound,g);	//TE mode
        tmeth.wait();	//multithread
    }
}

void simulator::low_TM_single_thread()
{
    for(uword i = 0; i < n; ++i)
    {
        g.iterate_region(reg,lapf);
        g.iterate_region(reg,relaxf);
        if(i % normalization_period == 0)
            normalize(1.0);
    }
	//JUST TESTING SOMETHING 12/27/2010
	//2,2 to 31,31
	write_file(g, "before.dat", grid::text, reg);
	grid::region nested;
	for(unsigned x = 2; x < 32; ++x)
		nested.push_back( grid::span(x,2,29) );
	functions::make_rand randomizer(29*29);
	g.iterate_region(nested,randomizer);
	settings::dt = 0.24;
	for(uword i = 0; i < n; ++i)
	{
		g.iterate_region(nested,lapf);
        g.iterate_region(nested,relaxf);
	}
	write_file(g, "after.dat", grid::text, reg);
	//END TESTING CODE 12/27/2010
}

void high_mode::subtract(grid::region & reg)
{
    grid & newm = modes.back();
    for(unsigned i = 0; i < modes.size()-1; ++i)
    {
        grid & oldm = modes[i];
        functions::inner_product ip( oldm );
        newm.iterate_region(reg, ip);
        grid::data_type inner_prod = tmeth.get_sum(ip.value());
        functions::sub_mode sm(oldm, sosq_vec[i], inner_prod);
        newm.iterate_region(reg, sm); //actually performs the subtraction
        tmeth.wait();
    }
}

//Note: high_mode::update_sumofsq() is NOT thread safe: it is called by main() only.
void high_mode::update_sumofsq(grid & g, grid::region & reg)
{
    functions::sum< square<grid::data_type> > sosq;
    g.iterate_region(reg, sosq);
    sosq_vec.push_back(sosq.value());
}

void simulator::normalize(grid::data_type maximum, grid::data_type scale)
{
    grid::data_type total_min = 0.0;
    if(!bound.empty())  //if TE mode
    {
        functions::min minime;
        g.iterate_region(reg, minime);
        total_min = tmeth.get_min(minime.value());
    }
    functions::max maximus;
    g.iterate_region(reg, maximus);
    grid::data_type total_max = tmeth.get_max(maximus.value());

    if(total_max - total_min < maximum)
    {
        if(std::abs(total_min) > maximum)
        {
            functions::add adder((total_min > 0) ? -total_min : total_min);
            g.iterate_region(reg, adder);
        }
        functions::multiply mult(scale * maximum / (total_max - total_min));
        g.iterate_region(reg, mult);
    }
}

//Sets each point in 'bound' equal to the avg of its partners
void TE_mode::adjust_boundary(const grid::boundary & bound, grid & gref)
{
    grid & g = gref;
    tmethod.wait();
    for(uword i = 0; i < bound.size(); ++i)
    {
        grid::data_type avg(0.0);
        for(uword j = 0; j < bound[i].second.size(); ++j)
            avg += g( bound[i].second[j].first, bound[i].second[j].second );
        avg /= static_cast<grid::data_type>(bound[i].second.size());
        g(bound[i].first.first, bound[i].first.second) = avg;
    }
}




